#!/usr/bin/python2.6

import sys
#import pickle
import string
import os
import io
import numpy as np
from scipy import stats
import time
#import epiCaseConconf

def snp_list_create(data_raw):
	"""Create a list of snps"""
	snp_list_merge = []		
	for epis in range(1,len(data_raw)):
		snp_list_merge.append(data_raw[epis][0])
		snp_list_merge.append(data_raw[epis][7])
	snp_list = list(set(snp_list_merge))
	snp_list_size = len(snp_list)#to be used in bonferoni
	return snp_list,snp_list_size

def cascon_calc_proper(jobCC):#,snplist_select):#snplist_select is an output file

	#sum of genotypes in file 1
	sum_genot1 = float(sum(jobCC["GenoList1"]))

	#sum of genotypes in file 2
	sum_genot2 = float(sum(jobCC["GenoList2"]))

	#sum of both files
	sum_g1g2 = float(sum_genot1 + sum_genot2)


	#colum sums
	sum_column = [jobCC["GenoList1"][0]+jobCC["GenoList2"][0],jobCC["GenoList1"][1]+jobCC["GenoList2"][1],jobCC["GenoList1"][2]+jobCC["GenoList2"][2]]

	#expected genotype frequences gent1. NB can be zero
	exp_genot1 = [(sum_genot1*sum_column[0])/sum_g1g2,(sum_genot1*sum_column[1])/sum_g1g2,(sum_genot1*sum_column[2])/sum_g1g2]

	#chisq-values for gent1
	if exp_genot1[0]>0:
		chig11=np.power((jobCC["GenoList1"][0] - exp_genot1[0]),2)/exp_genot1[0]
	else:
		chig11=0
	if exp_genot1[1]>0:
		chig12=np.power((jobCC["GenoList1"][1] - exp_genot1[1]),2)/exp_genot1[1]
	else:
		chig12=0
	if exp_genot1[2]>0:
		chig13=np.power((jobCC["GenoList1"][2] - exp_genot1[2]),2)/exp_genot1[2]
	else:
		chig13=0
	chig1=chig11+chig12+chig13

	#expected genotype frequences gent2
	exp_genot2 = [(sum_genot2*sum_column[0])/sum_g1g2,(sum_genot2*sum_column[1])/sum_g1g2,(sum_genot2*sum_column[2])/sum_g1g2]

	#chis-values for gent2
	if exp_genot2[0]>0:
		chig21=np.power((jobCC["GenoList2"][0] - exp_genot2[0]),2)/exp_genot2[0]
	else:
		chig21=0
	if exp_genot2[1]>0:
		chig22=np.power((jobCC["GenoList2"][1] - exp_genot2[1]),2)/exp_genot2[1]
	else:
		chig22=0
	if exp_genot2[2]>0:
		chig23=np.power((jobCC["GenoList2"][2] - exp_genot2[2]),2)/exp_genot2[2]
	else:
		chig23=0
	chig2=chig21+chig22+chig23	

	chi2_sum=float(chig1+chig2)

	chi2_pval= float(stats.chisqprob(chi2_sum,2))#final
	
	chiresult=(chi2_sum,chi2_pval)

	return chiresult	



